home *** CD-ROM | disk | FTP | other *** search
/ Mac Magazin/MacEasy 21 / Mac Magazin and MacEasy Magazine CD - Issue 21.iso / Wissenschaft & Technik / yorick12vr1-nofpu folder / include / spline.i < prev    next >
Text File  |  1995-07-26  |  10KB  |  308 lines
  1. /*
  2.    SPLINE.I
  3.    Cubic spline interpolator.
  4.  
  5.    $Id$
  6.  */
  7. /*    Copyright (c) 1994.  The Regents of the University of California.
  8.                     All rights reserved.  */
  9.  
  10. func spline(dydx, y, x, xp, dydx1=, dydx0=)
  11. /* DOCUMENT dydx= spline(y, x)
  12.        -or-   yp= spline(dydx, y, x, xp)
  13.        -or-   yp= spline(y, x, xp)
  14.      computes the cubic spline curve passing through the points (X, Y).
  15.  
  16.      With two arguments, Y and X, spline returns the derivatives DYDX at
  17.      the points, an array of the same length as X and Y.  The DYDX values
  18.      are chosen so that the piecewise cubic function returned by the four
  19.      argument call will have a continuous second derivative.
  20.  
  21.      The X array must be strictly monotonic; it may either increase or
  22.      decrease.
  23.  
  24.      The values Y and the derivatives DYDX uniquely determine a piecewise
  25.      cubic function, whose value is returned in the four argument form.
  26.      In this form, spline is analogous to the piecewise linear interpolator
  27.      interp; usually you will regard it as a continuous function of its
  28.      fourth argument, XP.  The first argument, DYDX, will normally have
  29.      been computed by a previous call to the two argument spline function.
  30.      However, this need not be the case; another DYDX will generate a
  31.      piecewise cubic function with continuous first derivative, but a
  32.      discontinuous second derivative.  For XP outside the extreme values
  33.      of X, spline is linear (if DYDX1 or DYDX0 keywords were specified,
  34.      the function will NOT have continuous second derivative at the
  35.      endpoint).
  36.  
  37.      The XP array may have any dimensionality; the result YP will have
  38.      the same dimensions as XP.
  39.  
  40.      If you only want the spline evaluated at a single set of XP, use the
  41.      three argument form.  This is equivalent to:
  42.           yp= spline(spline(y,x), y, x, xp)
  43.  
  44.      The keywords DYDX1 and DYDX0 can be used to set the values of the
  45.      returned DYDX(1) and DYDX(0) -- the first and last values of the
  46.      slope, respectively.  If either is not specified or nil, the slope at
  47.      that end will be chosen so that the second derivative is zero there.
  48.  
  49.      The function tspline (tensioned spline) gives an interpolation
  50.      function which lies between spline and interp, at the cost of
  51.      requiring you to specify another parameter (the tension).
  52.  
  53.    SEE ALSO: interp, tspline
  54.  */
  55. {
  56.   if (is_void(x)) {
  57.     /* spline(y,x) form */
  58.     x= y;
  59.     y= dydx;
  60.     dx= x(dif);
  61.     dy= y(dif);
  62.     diag= (2./dx)(pcen);
  63.     if (numberof(x)>2) diag(2:-1)*= 2.;
  64.     rhs= (3.*dy/(dx*dx))(pcen);
  65.     if (numberof(x)>2) rhs(2:-1)*= 2.;
  66.     dx= 1./dx;
  67.     if (is_void(dydx1)) {
  68.       if (is_void(dydx0)) {
  69.     return TDsolve(dx, diag, dx, rhs);  /* simple natural spline */
  70.       } else {
  71.     dx= dx(1:-1);
  72.     rhs= rhs(1:-1);
  73.     rhs(0)-= dydx0/dx(0);
  74.     return grow(TDsolve(dx, diag(1:-1), dx, rhs), dydx0);
  75.       }
  76.     } else {
  77.       dydx1= double(dydx1);
  78.       if (is_void(dydx0)) {
  79.     dx= dx(2:0);
  80.     rhs= rhs(2:0);
  81.     rhs(1)-= dydx1/dx(1);
  82.     return grow(dydx1, TDsolve(dx, diag(2:0), dx, rhs));
  83.       } else {
  84.     if (numberof(x)==2) return double([dydx1, dydx0]);
  85.     dx= dx(2:-1);
  86.     rhs= rhs(2:-1);
  87.     rhs(1)-= dydx1/dx(1);
  88.     rhs(0)-= dydx0/dx(0);
  89.     return grow(dydx1, TDsolve(dx, diag(2:-1), dx, rhs), dydx0);
  90.       }
  91.     }
  92.   }
  93.  
  94.   if (is_void(xp)) {
  95.     /* spline(y,x,xp) form */
  96.     xp= x;
  97.     x= y;
  98.     y= dydx;
  99.     dydx= spline(y,x,dydx1=dydx1,dydx0=dydx0);
  100.   }
  101.  
  102.   /* spline(dydx,y,x,xp) form */
  103.   l= digitize(xp, x); /* index of lower boundary of interval containing xp */
  104.   u= l+1;
  105.  
  106.   /* extend x, y, dydx so that l and u can be used as index lists */
  107.   dx= x(0)-x(1);
  108.   x= grow(x(1)-dx, x, x(0)+dx);
  109.   y= grow(y(1)-dydx(1)*dx, y, y(0)+dydx(0)*dx);
  110.   dydx= grow(dydx(1), dydx, dydx(0));
  111.  
  112.   xl= x(l);
  113.   dx= double(x(u)-xl);
  114.   yl= y(l);
  115.   dy= y(u)-yl;
  116.   dl= dydx(l);
  117.   du= dydx(u);
  118.   dydx= dy/dx;
  119.   return poly(xp-xl, yl, dl, (3.*dydx-du-2.*dl)/dx, (du+dl-2.*dydx)/(dx*dx));
  120. }
  121.  
  122. func tspline(tension, d2ydx2, y, x, xp, dydx1=, dydx0=)
  123. /* DOCUMENT d2ydx2= tspline(tension, y, x)
  124.        -or-     yp= tspline(tension, d2ydx2, y, x, xp)
  125.        -or-     yp= tspline(tension, y, x, xp)
  126.      computes a tensioned spline curve passing through the points (X, Y).
  127.  
  128.      The first argument, TENSION, is a positive number which determines
  129.      the "tension" in the spline.  In a cubic spline, the second derivative
  130.      of the spline function varies linearly between the points X.  In the
  131.      tensioned spline, the curvature is concentrated near the points X,
  132.      falling off at a rate proportional to the tension.  Between the points
  133.      of X, the function varies as:
  134.            y= C1*exp(k*x) + C2*exp(-k*x) + C3*x + C4
  135.      The parameter k is proportional to the TENSION; for k->0, the function
  136.      reduces to the cubic spline (a piecewise cubic function), while for
  137.      k->infinity, the function reduces to the piecewise linear function
  138.      connecting the points.  The TENSION argument may either be a scalar
  139.      value, in which case, k will be TENSION*(numberof(X)-1)/(max(X)-min(X))
  140.      in every interval of X, or TENSION may be an array of length one less
  141.      than the length of X, in which case the parameter k will be
  142.      abs(TENSION/X(dif)), possibly varying from one interval to the next.
  143.      You can use a variable tension to flatten "bumps" in one interval
  144.      without affecting nearby intervals.  Internally, tspline forces
  145.      k*X(dif) to lie between 0.01 and 100.0 in every interval, independent
  146.      of the value of TENSION.  Typically, the most dramatic variation
  147.      occurs between TENSION of 1.0 and 10.0.
  148.  
  149.      With three arguments, Y and X, spline returns the derivatives D2YDX2 at
  150.      the points, an array of the same length as X and Y.  The D2YDX2 values
  151.      are chosen so that the tensioned spline function returned by the five
  152.      argument call will have a continuous first derivative.
  153.  
  154.      The X array must be strictly monotonic; it may either increase or
  155.      decrease.
  156.  
  157.      The values Y and the derivatives D2YDX2 uniquely determine a tensioned
  158.      spline function, whose value is returned in the five argument form.
  159.      In this form, tspline is analogous to the piecewise linear interpolator
  160.      interp; usually you will regard it as a continuous function of its
  161.      fifth (or fourth) argument, XP.
  162.  
  163.      The XP array may have any dimensionality; the result YP will have
  164.      the same dimensions as XP.
  165.  
  166.      The D2YDX2 argument will normally have been computed by a previous call
  167.      to the three argument tspline function.  If you will be computing the
  168.      values of the spline function for many sets of XP, use this five
  169.      argument form.
  170.  
  171.      If you only want the tspline evaluated at a single set of XP, use the
  172.      four argument form.  This is equivalent to:
  173.           yp= tspline(tension, tspline(tension,y,x), y, x, xp)
  174.  
  175.      The keywords DYDX1 and DYDX0 can be used to set the values of the
  176.      returned DYDX(1) and DYDX(0) -- the first and last values of the
  177.      slope, respectively.  If either is not specified or nil, the slope at
  178.      that end will be chosen so that the second derivative is zero there.
  179.  
  180.      The function tspline (tensioned spline) gives an interpolation
  181.      function which lies between spline and interp, at the cost of
  182.      requiring you to specify another parameter (the tension).
  183.  
  184.    SEE ALSO: interp, tspline
  185.  */
  186. {
  187.   if (is_void(x)) {
  188.     /* tspline(tension, y,x) form */
  189.     x= double(y);
  190.     y= d2ydx2;
  191.     dx= x(dif);
  192.     dy= y(dif);
  193.     if (numberof(tension)==numberof(dx)) k= tension/abs(dx);
  194.     else k= tension*numberof(dx)/(max(x)-min(x));
  195.     k= max(min(k, 100./abs(dx)), 0.01/abs(dx));
  196.     kdx= k*dx;
  197.     skdx= sinh(kdx);
  198.     diag= (cosh(kdx)/skdx-1./kdx)/k;
  199.     diag= diag(pcen);
  200.     if (numberof(x)>2) diag(2:-1)*= 2.;
  201.     offd= (1./kdx-1./skdx)/k;
  202.     ddydx= (dy/dx)(dif);
  203.     if (is_void(dydx1)) {
  204.       if (is_void(dydx0)) {
  205.     if (numberof(x)==2) return [0., 0.];
  206.     diag= diag(2:-1);
  207.     offd= offd(2:-1);
  208.     return grow(0., TDsolve(offd, diag, offd, ddydx), 0.);
  209.       } else {
  210.     dydx0-= dy(0)/dx(0);
  211.     if (numberof(x)==2) return [0., dydx0/diag(0)];
  212.     diag= diag(2:0);
  213.     offd= offd(2:0);
  214.     ddydx= grow(ddydx, dydx0);
  215.     return grow(0., TDsolve(offd, diag, offd, ddydx));
  216.       }
  217.     } else {
  218.       dydx1= dy(1)/dx(1) - dydx1;
  219.       if (is_void(dydx0)) {
  220.     if (numberof(x)==2) return [dydx1/diag(1), 0.];
  221.     diag= diag(1:-1);
  222.     offd= offd(1:-1);
  223.     ddydx= grow(dydx1, ddydx);
  224.     return grow(TDsolve(offd, diag, offd, ddydx), 0.);
  225.       } else {
  226.     dydx0-= dy(0)/dx(0);
  227.     if (numberof(x)==2) return [dydx1/diag(1), dydx0/diag(0)];
  228.     ddydx= grow(dydx1, ddydx, dydx0);
  229.     return TDsolve(offd, diag, offd, ddydx);
  230.       }
  231.     }
  232.   }
  233.  
  234.   if (is_void(xp)) {
  235.     /* tspline(tension, y,x,xp) form */
  236.     xp= x;
  237.     x= y;
  238.     y= d2ydx2;
  239.     d2ydx2= tspline(tension, y,x,dydx1=dydx1,dydx0=dydx0);
  240.   }
  241.  
  242.   /* tspline(tension, d2ydx2,y,x,xp) form */
  243.   l= digitize(xp, x); /* index of lower boundary of interval containing xp */
  244.   u= l+1;
  245.  
  246.   /* extend x so that l and u can be used as index lists --
  247.      be careful not to make new intervals larger than necessary */
  248.   n= numberof(x)-1;  /* number of original intervals */
  249.   dxavg= (x(0)-x(1))/n;
  250.   if (dxavg>0.) {
  251.     dx0= max(max(xp)-x(0), dxavg);
  252.     dx1= max(x(1)-min(xp), dxavg);
  253.   } else {
  254.     dx0= min(min(xp)-x(0), dxavg);
  255.     dx1= min(x(1)-max(xp), dxavg);
  256.   }
  257.   x= grow(x(1)-dx1, x, x(0)+dx0);
  258.  
  259.   /* compute k so that sinh(k*dx) is safe to compute */
  260.   dx= x(dif);
  261.   if (numberof(tension)==n) {
  262.     k= grow(0., tension, 0.)/abs(dx);
  263.   } else {
  264.     k= tension/abs(dxavg);
  265.   }
  266.   k= max(min(k, 100./abs(dx)), 0.01/abs(dx));
  267.  
  268.   /* extend y carefully so that linear extrapolation happens automatically */
  269.   k1= k(2);
  270.   k0= k(-1);
  271.   dydx1= (y(2)-y(1))/dx1;
  272.   kdx= k1*dx1;
  273.   d2u= d2ydx2(2);
  274.   d2l= d2ydx2(1);
  275.   dydx1+= ((d2u-d2l*cosh(kdx))/sinh(kdx) - (d2u-d2l)/kdx)/k1;
  276.   dydx0= (y(0)-y(-1))/dx0;
  277.   kdx= k0*dx0;
  278.   d2u= d2ydx2(0);
  279.   d2l= d2ydx2(-1);
  280.   dydx0+= ((d2u*cosh(kdx)-d2l)/sinh(kdx) - (d2u-d2l)/kdx)/k0;
  281.  
  282.   y= grow(y(1)-dydx1*dx1, y, y(0)+dydx0*dx0);
  283.   d2ydx2= grow(0., d2ydx2, 0.);
  284.  
  285.   /* begin interpolation */
  286.   xu= x(u);
  287.   xl= x(l);
  288.   dx= xu-xl;
  289.   dxl= xp-xl;
  290.   yl= y(l);
  291.   dydx= (y(u)-yl)/dx;
  292.  
  293.   km2= 1./(k*k);
  294.   km2(1)= 0.;
  295.   km2(0)= 0.;
  296.   km2= km2(l);
  297.   k= k(l);
  298.   skdx= sinh(k*dx);
  299.  
  300.   d2u= d2ydx2(u);
  301.   d2l= d2ydx2(l);
  302.   d3= km2*(d2u-d2l)/dx;
  303.  
  304.   d2ydx2= d2u*sinh(k*dxl)/skdx + d2l*(sinh(k*(xu-xp))/skdx-1.);
  305.  
  306.   return yl + km2*d2ydx2 + (dydx-d3)*dxl;
  307. }
  308.